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' In this paper, a lattice Boltzmann (LB) model is presented for axisymmetric multiphase flows. 

' Source terms are added to a two-dimensional standard lattice Boltzmann equation (LBE) for multi- 

, phase flows such that the emergent dynamics can be transformed into the axisymmetric cylindrical 

coordinate system. The source terms are temporally and spatially dependent and represent the 
' axisymmetric contribution of the order parameter of fluid phases and inertial, viscous and surface 

tension forces. A model which is effectively explicit and second order is obtained. This is achieved by 
taking into account the discrete lattice effects in the Chapman-Enskog multiscale analysis, so that 
QQ ' the macroscopic axisymmetric mass and momentum equations for multiphase flows are recovered 

self-consistently. The model is extended to incorporate reduced compressibility effects. Axisym- 
metric equilibrium drop formation and oscillations, breakup and formation of satellite droplets from 
viscous liquid cylindrical jets through Rayleigh capillary instability and drop collisions are presented. 
C ' Comparisons of the computed results with available data show satisfactory agreement. 

>»: 

. PACS numbers: 47.11. 47.55.Kf, 05.20.Dd, 47.55. Dz, 47.20.Ma 
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Fluid flow with interfaces and free surfaces is common in nature and in many engineering applications. Such 
interfacial flows which typically involve multiple scales remain a formidable non-linear problem rich in physics and 
continue to pose challenges to experimentalists and theoreticians alike . Numerical simulation of multiphase flows 
is challenging as the shape and location of the interfaces must be computed in conjunction with the solution of the 
flow field |3|| . Computational methods based on the lattice Boltzmann equation (LBE) for simulating complex 
emergent physical phenomena have attracted much attention in recent years ]3 1^ . The LBE simulates multiphase 
' flows by incorporating interfacial physics at scales smaller than macroscopic scales. Phase segregation and interfacial 
7-H fluid dynamics can be simulated by incorporating inter-particle potentials 6, 7J, concepts based on free energy 

■ or kinetic theory of dense fluids 10; I L, id . 

■ The formulation of the standard LBE is based on the Cartesian coordinate system and does not take into account 
jn , axial symmetry that may exist. Numerous multiphase flow situations exist where the fluid dynamics can be ap- 

■ proximated as axisymmetric [H Il3||. Examples include head-on collision of drops, normal drop impingement on solid 
^ : surfaces and Rayleigh instability of cylindrical liquid columns. Currently, full three-dimensional (3D) calculations 

. ^ have to be carried out for problems which may be approximated as axisymmetric 0,0,^^. In 3D computations, 
^ \ computational considerations restrict the numerical resolution that may be employed and the physics may not be 
well resolved. For example, in breakup of drops into satellite droplets the size of the droplets may be such that the 
3D grids may not resolve them. To improve the computational efficiency of the LBE for axisymmetric multiphase 
flows, we propose an axisymmetric LB model in this paper. The approach consists of adding source terms to the 
.„ two-dimensional (2D) Cartesian LBE model based on the kinetic theory of dense fluids for multiphase flows [lol[ll| . 
, This approach is similar in spirit to the idea proposed in [17| to solve single-phase axisymmetric flows. However, 
^ ' multiphase flow problems involve additional complexity as a result of interfacial physics involved, i.e. the surface 
J tension forces and the need to track the interfaces. In this case, the accuracy of the numerical discretization of the 
source terms representing interfacial physics also becomes an important consideration. 

This paper is organized as follows. In Section m the axisymmetric LBE multiphase model is described. Then, in 
Section [nil its extension to simulate axisymmetric multiphase flows with reduced compressibility effects is described. 
The computational methodology adopted is also discussed in this section. In Section llVI the axisymmetric model is 
applied to benchmark problems to evaluate its accuracy. Finally, the paper closes with summary in Section Ivi 
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II. AXISYMMETRIC LBE MULTIPHASE FLOW MODEL 



To simulate axisymmetric multiphase flows, axisymmetric contributions of the order parameter, and inertial, viscous 
and surface tension forces may be introduced to the standard 2D LBE. The source terms, which will be shown to be 
spatially and temporally dependent, are determined by performing a Chapman-Enskog multiscale analysis in such a 
way that the macroscopic mass and momentum equations for multiphase flows are recovered self-consistently. The 
introduction of source terms makes it necessary to calculate additional spatial gradients when compared to those in 
the standard LBE. While this approach is developed for a specific LBE multiphase flow model based on kinetic theory 
of dense fluids 0,0], it can be readily extended to other LBE multiphase flow models. 

The governing continuum equations of isothermal multiphase flow [T8[ll9j| in the cylindrical coordinate system when 
the axisymmetric assumption is employed are 



dtp + -dr (prur) + {pu^) = 0, 
r 



(1) 



p (dtUr + UrdrUr + U^dzUr) = -drP + Fs.r + Fext,r + -dr {rlirr) + dz (Ilrz) , 

r 



(2) 



P (dtUz + UrdrUz + UzOzUz) 



'dzP + Fs.z + Fe, 



1 



dr (rllzr) + dz (Uzz) ., 



(3) 



where p is the density and Ur and Uz are the radial and axial components of velocity. These equations are derived 
from kinetic theory that incorporates intermolecular interactions forces which are modeled as a function of density 
following the work of van der Waals |23| . The exclusion volume effect of Enskog j^J] is also incorporated to account 
for increase in collision probability due to the increase in the density of non-ideal fluids. These features naturally give 
rise to surface tension and phase segregation effects. The other variables which appear in the above equations will 
now be described. Hrr, ^rz, ^zz are the components of the viscous stress tensor and are given by 



Urz = lizr = P (dzUr + drUz) ■ 

Uzz = 2pdzUz, 



(4) 
(5) 
(6) 



where p is the dynamic viscosity. F^ r and F^ z are the axial and radial components respectively of the surface tension 
force, which are given by |19| 



Fs^r — Updr 

Fs^z = updz 



dr{rdrp) + dz{dzp) 



-dr{rdrp) + dz{dzp) 



(7) 
(8) 



where k controls the strength of the surface tension force. This parameter is related to the surface tension of the fluid, 
CT, through the density gradient across the interface by the equation 



dp 
dn 



dn. 



(9) 



Thus, the surface tension is a function of both the parameter k, and the density profile across the interface. The 
terms Fext,r and F^xt.z in Eqs. Q and (j2Jl respectively are the radial and axial components of external forces such as 
gravity. 

The pressure, P, is related to density through the Carnahan-Starling-van der Waals equation of state (EOS) 



P = pRT 



1 + 7 + 7^-7^ 



ap 



(10) 



where 7 = hp/ A. The parameter a is related to the intermolecular pair- wise potential and h to the effective diameter of 
the molecule, d, and the mass of a single molecule, m, by 6 = 2ii(P /ira. i? is a gas constant and T is the temperature. 
The Carnahan-Starling EOS has a supernodal P — 1/ p — T curve, i.e., dP/dp < 0, for certain range of values of p, 
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FIG. 1: Schematic of arrangement of coordinate system in axisymmetric multiphase flow ((r, z) and {y, x) coordinate directions 
are shown). 

when the state fluid temperature is below its critical value. This unstable part of the curve is the driving mechanism 
responsible for keeping the phases of fluids segregated and for maintaining a self-generated sharp interface. 

We now modify the standard LBE in such a way that it effectively yields the axisymmetric multiphase flow equations, 
Eqs. (0- H10|l . in a self-consistent way. To facilitate this, we employ the following coordinate transformation, 
illustrated in Fig. Q which allows the governing equations to be represented in a Cartesian-like coordinate system, 
i.e. {x,y): 

{r,z)^{y,x), (11) 



Assuming summation convention for repeated subscript indices, Eqs. ((T|l- 

y 



may be transformed to 



dtp + dk [puk] 



(12) 



(13) 



where 



p {dtUi + UkdkUi) = -diP -I- Fs^i + Fext,t + dk [p [dkUi + diUk)] + Fax,i, 



(14) 



(15) 



and k £ {x, y}. The right hand side (RHS) in Eq. (|13|1 . —puy/y, is the additional term in the continuity equation 
that arises from axisymmetry. The corresponding term for the momentum equation, Eq. (|14|) . is 



F, 



ax.i ^ ^ [dyUi + diUy] -t- npdi \ ~^'^vP 



(16) 



To recover Eqs. H13(l and 1)14(1 . we introduce two additional source terms, and S^, to the standard 2D Cartesian 
LBE which has U,a as its collision term and a source term for the internal and external forces, Sa- These unknown 
additional terms, representing the axisymmetric mass and momentum contributions respectively, are to be determined 
so that the macroscopic behavior of the proposed LBE corresponds to axisymmetric multiphase flow. Thus, we propose 
the following LBE 

fa{x + ea5t,t^- St) - fa{x,t) = '^[^a\{x .t) + ^a\{x+e^5t ,t+St)\ + 



[Sa\{x,t) + Sa\(x+ec,St.t+St)\ + 
Sa\ix,t) + Sa\{x+eo,St,t+St) + 
Sa\{x.t) + Sa\{x+ec,5t,t+St) ^t, 



(17) 
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where fa is the discrete single-particle distribution function, corresponding to the particle velocity, e^, where a is 
the velocity direction. The Cartesian component of the particle velocity, c, is given by c = S^/St, where S.j; is the 
lattice spacing and 5t is the time step corresponding to the two-dimensional, nine-velocity model(D2Q9) ^24] shown 
in Fig. n Here, the collision term is given by the BGK approximation 

= r^^, (18) 

T St 

where A is the relaxation time due to collisions, 5t is the time step and /^'^ is the truncated discrete form of the 
Maxwellian 

pg^pg.Af/ N J Ca-M {eg- uf 1 U- U \ 

/a -/o + 2(i?r)2 2 RT j ' ^ 

where R is the gas constant, T is the temperature and Wa is the weighting coefficients in the Gauss-Hermite quadrature 
to represent the kinetic moment integrals of the distribution functions exactly 26] . For isothermal flows, the factor 
RT is related to the particle speed c as RT = l/Sc^. The term in Eq. ifTTI) 

= (^"^-"^-)(J}+fe.u) _^e,,M(^^^) (20) 

pRl 

represents the effect of internal and external forcing terms on the change in the distribution function. The internal 
force term gives rise to surface tension and phase segregation effects which are given by 

Fj=-dj^ + Fs,j, (21) 

where the function ip = P — pRT is the non- ideal part of the equation of state given in Eq. (|10|l . The first two terms 
on the RHS of Eq. (|17|l corresponds to those presented by He et al. (1998). As mentioned above, the last two terms, 
and , in this equation is to be selected such that its behavior in the continuum limit would simulate the influence 
of the non-Cartesian-like terms in Eqs. (|13|l and H14|) in a self-consistent way. Since the zeroth kinetic moment of 
the term fa'^'^^{p,0) is involved in the derivation of the macroscopic mass conservation equation from the LBE, the 
source term in Eq. (|17|l is proposed to be equal to f^'^\p,{)) multiplied by an unknown m and normalized by 
the density p. The other source term 5',^ is proposed analogous to the source term in Ea. (|2l)|) . Thus, we propose 

p 

q" _ (ggj - req,M 



pRT 



rj''"iP.u)- (23) 



Here the unknowns, m and , in the above two equations can be determined through Chapman-Enskog analysis as 
will be shown later. It must be stressed that all terms, including the collision term, on the RHS are discretized by 
the application of the trapezoidal rule, since it has been arg ued that at least a second-order treatment of the source 
terms is necessary for simulation of multiphase flow |lOl lll| . The macroscopic fields are given by 

P = E/- (24) 

a 

pUi = E/aCQi. (25) 

a 

In this model, the order parameter is the density, p, which distinguishes the different phases in the flow. 

Equation (|17|l is implicit in time. To remove implicitness in this equation we introduce a transformation following 
the procedure described by He and others 0, , whereby 



2 " 2 
in Eq. H17|) . so that we obtain 

- - T 

fa{x + eaSt,t + 6t) - fa{x,t) = ila\(x,t) H — TTTT 

T -t- i/Z 



fa=fa-l:na-USc.+Sa + Sa]St (26) 



Sa + S„ 



\{x,t)"t 



5u (27) 



5 



where 

= -l^. (28) 

Thus, fa is the transformed distribution function that removes impUcitness in the proposed LBE, Eq. (|17|l . which 
describes the evolution of the fg distribution function. The foUowing constraints on the equihbrium distribution and 
the various source terms [2^ are imposed from their definition: 

a ct a 

fa'^SaieajCak = P{RTY {u^Sjk + UjSki + UkSij) , (29) 

a 

^ = 0, ^ S'aeai = , ^ SaCaiCaj = [Fi + Fext,i)Uj + (Fj + Fext.j)Ui, (30) 
a a a 

^ 5^ = m , ^ S'^ecj = 0, ^ S^eateaj = m i?r5y, (31) 

a a a 

^Sa^O, '^a^"* ^ ' ^c^eazeaj = (F, + Fj Ui). (32) 

a a a 

Then the following relationships are obtained between the transformed distribution function and the macroscopic 
fields, which also include the curvature effects resulting from axial symmetry: 

P = Yfa + ^rn'St, (33) 

a 

pu^ = ^/ae„, + i(F, +Fe,t,, + F;>t. (34) 

a 

Now, to establish the unknowns m and F^ in the above formulation, the Chapman-Enskog multiscale analysis is 
performed 21]. Introducing the expansions |30| 



fa{x + ea6t,t + 6t) = J2^tJgix,t), (35) 

0=0 

A„ = dt„ +eakdk, (36) 

OC 



oo 



dt = Ee"at„, (38) 



0=0 



where e = St in Eq. (|27|l and using Eq. (|26|l to transform fa back to /„, the following equations are obtained in the 
consecutive order of the parameter e: 

0(e°) : fi"^ = fl\ (39) 
0(ei) : Ao/^"^ = -^/^'^ + 5„ + 5; + 5;:, (40) 

0(.^):a,j("'+Ao/« = -^/P. (41) 
Now, invoking the Chapman-Enskog ansatz 



(42) 
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and performing J2ai') '-'^ ^is. H40() and H41|l . we obtain 

dtoP + dkipuk) = m, (43) 
dt,P = 0, (44) 

respectively. Combining the first- and second- order results given by Eqs. H43|) and 144|l and considering dt = dto +£f^ti j 
we get 

dtp + dkipuk) ^ m . (45) 
Comparing this equation and Eq. the unknown m is obtained as 

m' = (46) 

y 

This is the axisynimetric contribution to the Cartesian form of the equation for the order parameter, i.e., density 
characterizing the different phases of the flow. Taking the first kinetic moment, J2a^ai{'), of Eqs. (|40|l and 141|l . 
respectively, we get 

dtoipui) + dkipu^Uk) = -d,{pRT) + F, + F,,t., + F- , (47) 

dtApui) + dkUlf - 0, (48) 

where 

n.?=E/^'^^-^-^- (49) 

a 

Employing the expression for /q^-* from Eq. H40() in Eq. H49() . together with the summational constraints given above, 
and neglecting terms of the order 0{Ma^) or higher, we get 

n|j^ = -TRTp{djU, -t- ^^UJ). (50) 

Equation then simplifies to 

dtt (pui) = dj {TRTp{djUi + diUj)) . (51) 

Combining Eqs. 147|l and l|51|l . we get 

dtipu^) + dk{pu,Uk) = -d^ipRT) + F, + F^^.^ + F- + d, {rStRTpidjU, + d,u/)) , (52) 
or substituting for Fi from Eq. H21|l . we obtain 

dtipu,) + dk{pu,uk) = -d^P + F,,, + F,,t^, + F'i + d, {T5tRTp{djU^ + ^^Uj)) . (53) 
Using Eqs. (|45() and H53() . this can be simplified to 

p{dtu, + ukdku,)-^^ - -a,P + F,,, + J^e^M + + 

y 

dj {T5tRTp{djU, + d,uj)) . (54) 
Comparing Eqs. (|14|) and H54|) . we obtain the other unknown F.^ where 

f'; = Fa,, ~ ^ = ^ [dyu, + d,uy] + npd, Qa.p) - (55) 

This is the axisymmetric contribution to the Cartesian form of the equation for the momentum, where the first, 
second and the third terms on the RHS correspond to the viscous, surface tension and inertial force contributions, 
respectively. The dynamic viscosity is related to the relaxation time for collisions hy p = prStRT = pXc^, where 
= l/3c^. The set of equations corresponding to the axisymmetric LBE multiphase flow model is given by Eqs. 
(|?7jl and together with Eqs. (jSOJ, l|23 and (j2Sll, 1231 and (EH, and (gHl and (jSHll- In general, this muhiphase 
model and that proposed by He and others 10] face difficulties for fluids far from the critical point and/or in the 
presence of external forces. This difficulty is related to the calculation of the intermolecular force in Ea. (|21|l . involving 
the computation of djijj which can become quite large across interfaces. Unless this term is accurately computed, 
the model may become unstable because of numerical errors |l4l IbH . Hence, an improved treatment of this term is 
necessary. This will now be described. 
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III. AXISYMMETRIC LBE MULTIPHASE FLOW MODEL WITH REDUCED COMPRESSIBILITY 

EFFECTS 

He and co-workers have proposed that through a suitable transformation of the distribution function, fa , which 
involves invoking the incompressibility condition of the fluid, and employing a new distribution function for capturing 
the interface, the difficulty with handling the intermolecular force term, djijj, can be reduced. We apply this idea 
to the axisymmetric model developed in the previous section. We replace the distribution function fa by another 
distribution function through the transformation [T]| 



P 



(56) 



The effect of this transformation will be discussed in greater detail below. By considering the fluid to be incompressible, 
i.e. 



at 



and using the transformation Eqs. H56|) and ^2<6\ . Eq. (|27|) is replaced by 



ga{x + eQ(5t,t + 5t) - ga{x,t) = ^ga\{x,-: 



where 



and 



r + 1/2 



Sqa + 



ga 



S' 



ga 



^qa — 



gl' = fa'RT + ^{p) 



9a - 97 
T+1/2' 

f^''^'ip,0) 



(57) 



(58) 



(59) 



(60) 



The corresponding source terms become 

Sga — {Saj ~ Uj) X 



(p^p ,f^:^lM ( ra''^{p.u) /^^'^^(p,o) 1^ , 

{Fj + Fe^tj) OjVAP) 



(61) 



^ga — SaRT 



0) 



y J 



(62) 



Sqa =SaRT ^ [Cj - Uj)Fj 



ra'^^'{p,u) 



(63) 



The term in Eq. H61() is multiplied by the factor (/^*'*^(/0, u)/ p — /q*'*^(p, 0)//o). This factor, from the definition 
of the equilibrium distribution function, f^, in Eq. (|19|l is proportional to the Mach number and thus becomes 
smaller in the incompressible limit. Hence, it alleviates the difficulties associated with the calculation of the djij}, a 
major source of numerical instability with the original model 10] . Thus, Eqs. H58|l - H63l) are found to be numerically 
more stable compared to Eq. H27|l supplemented with Eqs. (|20() . H22I) and (|23l) . In this new framework, we still need 
to introduce an order parameter to capture interfaces. Here, we employ a function, 0, referred to henceforth as the 
index function, in place of the density, as the order parameter to distinguish the phases in the flow. 

The evolution equation of the distribution function whose emergent dynamics govern the index function has to 
be able to maintain phase segregation and mass conservation. To do this, we employ Eq. (|27() together with Eqs. 
(|20|l . (|22|) and H23|l by keeping the term involving djij: and m , while the rest of the terms may be dropped as they play 
no role in mass conservation. In addition, the density is replaced by the index function in these equations. Hence, 
the evolution of the distribution function for the index function is given by 



fa{x + eaSt,t + St) - fa{x,t) = ^lfa\{x,t} + 



1/2 



S 



fa 



S 



fa 



l(x,t)t 



(64) 
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where the collision and the source terms are given by 



f — S. f 



p p \ y 

The hydrodynamic variables such as pressure and fluid velocity can be obtained by taking appropriate kinetic moments 
of the distribution function ga, i.e. 

P = + (68) 

a 

pRTui = ^ gaCai + ^ {Fs^i + Fe^t,i) (5t + (5t. (69) 

a 

This follows from the definition of ga given in Eq. (|56|l and also includes curvature effects. The index function is 
obtained from the distribution function by taking the zeroth kinetic moment, i.e. 

'^ = E/" + ^-™'^*- (70) 
^ — ' 2 p 

The terms ra and are given in Eqs. H46|l and H55|l. respectively. The density is obtained from the index function 
through linear interpolation, i.e. 

p{4>)=PL + 4- — '^(ph-Pl), (71) 
<Ph - <Pl 

where pL and pu are the densities of the light and heavy fluids, respectively, and 0^ and (pn refer to the minimum 
and maximum values of the index function, respectively. These limits of the index function are determined from 
Maxwell's equal area construction applied to the function ip{4') + 4>RT. 

Thus, the axisymmetric LBE multiphase flow model with reduced compressibility effects corresponds to Eqs. (|58f) - 
(|71|l . The relaxation time for collisions is related to the viscosity of the fluid using the same expression as derived in 
the previous section. If the kinematic viscosity of the light fluid, v^, is different from that of the heavy fluid, vu, its 
value at any point in the fluid is obtained from the index function through linear interpolation, i.e. 

^('^)=^+/^(^-'^^)- (72) 
<Ph - <Pl 

It may be seen that the model requires the calculation of spatial gradients in Eqs. and (|66|l and of the 

Laplacian in Eq. (|15|l . Since maintaining accuracy as well as isotropy is important for the surface tension terms, they 
are calculated by employing a fourth-order flnite-difference scheme for the gradient and a second-order scheme for the 
Laplacian, given respectively by 



^^^ = ^Y. [^^(^ + ^-'^*) - ^(^ + '^^»^5t)] (^) + 0{5t), (73) 

^ a — 1 

and 



1 ^ 

V^tu = ^^^^w = ^Y^ [^{x + e„,<5t) - w{x)] -f 0[5l), (74) 

^ a=l 

for any function m. Notice that these discretizations are both based on the lattice based stencil, instead of the 
standard stencil based on the coordinate directions. In addition, in the application of this model, the implementation 
of boundary conditions plays an important role. In particular, along the axisymmetric line, i.e. y = 0, specular 
reflection boundary conditions are employed for the distribution functions. For the two-dimensional, nine velocity 
(D2Q9) model shown in the inset of Fig. [U we set /a = /4, /g = /§, /g = f^ and 52 = 54, 55 = 58 and g^ = 57 
for the distribution functions after the streaming step. For macroscopic conditions, along this line, Uy = dy{-) = 0, 
through which the singular source terms of type l/y(- ) in the model can be appropriately treated. On the other hand, 
boundary conditions along the other lines are similar to those for the standard LBE. 
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FIG. 2: Pressure difference across axisymmetric drops as a function of radius for different values of the surface tension parameter 
k; Comparison of computed results using the axisymmetric LBE model versus theoretical prediction based on the Laplace- Young 
relation. Quantities are in lattice units. 



IV. RESULTS AND DISCUSSION 



In the rest of this paper, unless otherwise specified, the results are presented in lattice units, i.e. the velocities are 
scaled by the particle velocity c, the distance by the minimum lattice spacing 6^ and time by c/6x. All other quantities 
are scaled as appropriate combinations of these basic units. First, the axisymmetric LBE multiphase flow models 
are applied to verify the well-known Laplace- Young relation for an axisymmetric drop. According to this relation, 
AP = 2a/Rd, where AP is the difference between the pressure inside and outside of a drop, a is the surface tension 
and Rd is the drop radius. For different choices of the surface tension parameter, k, the surface tension values are 
obtained from Eq. (jSJ by the replacing density in Eq. lO and (jSJ by the index function. To obtain the normal gradient 
used in Eq. a physical configuration consisting of a liquid and a gas layer is set up. Once equilibrium is reached, 
the density gradient may be computed and hence the surface tension. Having obtained the relationship between the 
surface tension cr, and the parameter k, axisymmetric drops of four different radii, Rd = 40, 50, 60 and 70, are set up 
in a domain discretized by 201 x 101 lattice sites. Periodic boundaries are considered in the x direction and an open 
boundary condition is considered along the boundary that is parallel to the axisymmetric boundary. By considering 
three different values of k, 0.05, 1.0 and 0.15, the pressure difference across the drops is determined. Figure [21 shows 
a comparison of the pressure difference across the interface of the drops computed using the axisymmetric model 
developed in Section UTTl and that predicted by the Laplace- Young relation. It is found that the computed results are 
in good agreement with the theoretical values, with a maximum relative error of about 3%. 

Another important test problem is that of an oscillating axisymmetric drop immersed in a gas. Since current 
versions of the LBE simulate a relatively viscous fluid, it is appropriate to compare the oscillation frequency with 
that of Miller and Scriven (1968) . In contrast to earlier analytical solutions on drop oscillations, this work considers 
viscous dissipation effects in the boundary layer at the interface. According to the frequency for the n*'' mode 
of oscillation for a drop is given by 

UJn = LJn ~ ^CtlVn^ + ^Q;^, (75) 

where a;„ is the angular response frequency, and w* is Lamb's natural resonance frequency expressed as 

2 _ n(n + l)(n-l)(n + 2) 
^""^ - i?3[np, + (n+l)p.] " ^''^ 



Rd is the equilibrium radius of the drop, a is the interfacial surface tension, and pi and pg are the densities of the two 
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FIG. 3: Configurations of an oscillating drop as a function of time; Rmin ~ 40, Rmax = 55, pg = 0.1, pi = 0.4, vi = Vg 
1.6667 X 10^^. Quantities are in lattice units. 



fluids. The parameter a is given by 



a — 



2iRd[npg + {n + l)pi] {pipi)l + {pgPg) 



(77) 



where /i; and pg are the dynamic viscosity of the two liquids. The subscripts g and I refer to the ambient gas and 
liquid phases, respectively. We consider the second mode of oscillation and analytical expressions for the time period 
are presented in Eq. (|75|l . 

The initial computational setup consists of a prolate spheroid of minimum {Rmin) and maximum (Rmax) radii of 40 
and 55, respectively, placed in the center of the domain discretized by 201 x 101 lattice sites. We consider the surface 
tension parameters: k — 0.2, and the density of the gas and the drop to be pg = 0.1 and pi = 0.4, respectively. The 
kinematic viscosity of both the gas and the drop are considered to be the same and given by Ug = ui ^ 1.6667 x 10'^. 
Figure O shows the configurations of an oscillating drop at different times computed using the standard axisymmetric 
model with these conditions. The drop changes from a prolate shape at i = 2000 to oblate shape at i = 16000. Such 
shape changes continue till the drop reaches its equilibrium spherical shape. Figure 0] shows the temporal evolution 
of the interface locations of the oscillating drop with the conditions above for two different surface tension parameter: 
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FIG. 4: Interface location of an oscillating drop as a function of time for two values of the surface tension parameter k; 
Rmin = 40, Rmax = 55, pg — 0.1, pi = 0.4, vi — Vg = 1.6667 X 10~^. Quantities are in lattice units. 



K = 0.02 and 0.08. It is expected that increasing the surface tension will reduce the time period of oscillations. The 
computed (Tlbe) and analytical (Tanai) time periods, where Tanai = 2^/102, when k = 0.02 are 29483 and 29448 
respectively. As k is increased to 0.08, Tlbe and Tanai become 14388 and 14313 respectively. It may be seen that 
the computed and analytical values agree well, the difference being less than 1%. Also, the time period decreases as 
K is increased, which is consistent with expectations. 

Consider next the effect of changing the drop size on the time period of oscillations. Figure [51 shows the interface 
locations of an oscillating drop as a function of time for the following two initial sizes: Rmin = 30 and Rmax — 45; 
Rmin — 40, Rmax — 55. Reducing the drop size reduces its time period. The computed time period of the larger 
drop is equal to 29483, while that for the smaller drop is 20118. Comparison of the computed time periods with the 
analytical solution shows that they agree within 1% for these cases. Next, consider three different kinematic viscosities 
of the hquid: i^i — 1.6667 x 10^^,3.3333 x 10~^ and 5.0 x 10^^. Figure shows the effect of drop viscosity on the 
temporal evolution of the interface locations of the drop. It is found that as the kinematic viscosity is increased the 
time period increases moderately which is consistent with the analytical solution. The computed time periods at these 
viscosities are 29483, 31030 and 32925, while the analytical values are 29448, 30597 and 31318, respectively, with a 
maximum error within 5.1%. 

The third test problem considered here is that of the break-up of a cylindrical liquid column into drops, a fascinating 
problem of long standing theoretical and practical interest. In a seminal work, Rayleigh (1878) showed through a 
linear stability analysis of an inviscid column of cylindrical liquid of radius Rc that the column will be unstable if the 
axisymmetric wavelength of any disturbance is longer than its circumference, i.e. the wave number k* = 2'KRc/Xd 
should be less than one. Later, the theoretical analysis was extended to more realistic conditions by including viscosity. 
In the last three decades, several experimental and numerical investigations have also been performed. To evaluate 
the axisymmetric LBE model, we study the Rayleigh capillary instability for different wavenumbers. Initial studies 
carried out with k* > 1 showed that the liquid does not break-up. We will now present results of cases with break-up. 
Consider a cylindrical liquid column of radius Rc = 45 subject to an axisymmetric co-sinusoidal wavelength \d — 320, 
i.e. k* — 0.88. To simulate the dynamics of instability for this wavenumber, we consider a domain discretized by 
321 X 151 lattice sites with pg = 0.1, pi = 0.4, Vg = vi = 6.6667 x 10~^ and k — 0.1. Since k* < 1, it is expected that 
the liquid column would eventually breakup. FigureElshows the configurations of the liquid column at different times. 
As time progresses, the imposed interfacial disturbances on the liquid column grow. At t — 28000, 46000 52000, the 
cross-section of the column becomes progressively thinner in the center, and by mass conservation, the ends becomes 
larger. At t = 60000, notice that a bead-type structure is formed at the ends and with a thin ligament between them. 
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FIG. 5: Interface location of an oscillating drop as a function of time for two drop sizes; Pg = 0.1, pi = 0.4, vi = Vg = 
1.6667 X IQ-^ K = 0.02. Quantities are in lattice units. 
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FIG. 6: Interface location of an oscillating drop as a function of time for different kinematic viscosities vi; Rmin = 40, Rmax = 55, 
Pg = 0.1, pi = 0.4, K = 0.02. Quantities are in lattice units. 
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FIG. 7: Configurations of a cylindrical liquid column at difi'erent times undergoing Rayleigh breakup and satellite droplet 
formation; k* = 0.88, Pg — 0.1, pi = 0.4, vi — Vg — 6.6667 x 10~^. Quantities are in lattice units. 



Such a structure has been observed in experiments m and in other numerical simulations |35l | . Eventually, the column 
breaks up forming a thin ligament in the middle, which then becomes a satellite droplet. 

Let us now increase the wavelength of the disturbance to = 600, keeping the physical parameters the same as 
before. We consider a domain represented by 601 x 151 lattice sites. Since, Rc — 45 as before, the wavenumber is 0.47. 
Figure 13 shows the temporal evolution of the configurations of the liquid column at this reduced wavenumber. The 
axisymmetric disturbance grows with time. Since the wavelength is longer, it can be noticed that the ligament that is 
formed during the Rayleigh instability is also longer. As a result, after the column breaks up, a larger satellite droplet 
is formed. To express the drop size distribution with wave numbers more quantitatively, we plot the non-dimensional 
size of the main and satellite drops, r* = R/Rc, as a function of wave number, k* in Fig. |3 It may be noted that 
Rayleigh's original analysis predicts only the onset of breakup and not the formation of satellite droplets. To predict 
analytically satellite droplet formation, it has been shown that at least a third-order perturbation analysis of the 
Navier-Stokes equations (NSE) is needed [s^. Computations based on direct solutions of the NSE also predict the 
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FIG. 8: Configurations of a cylindrical liquid column at diff'erent times undergoing Rayleigh breakup and satellite droplet 
formation; k* — 0.47, pg — 0.1, pi = 0.4, vi — Ug — 6.6667 x 10~^. Quantities are in lattice units. 



formation of the satellite droplets. 

To evaluate the drop size distribution computed using the axisymmetric LBE model, we consider the experimental 
data of Rutland and Jameson (1971), the experimental data and analytical solution based on a third-order perturbation 
analysis of the NSE by Lafrance (1975), a boundary integral solution of the NSE by Mansour and Lundgren (1990) and 
a finite element solution of the NSE by Ashgriz and Mashayek (1995). It can be seen in the figure that as long as the 
wavenumber is less than one, as expected there will be a satellite droplet formation. As the wavenumber is reduced, 
the sizes of both the main drop and satellite droplet increases. The rate of increase of the size of the satellite droplet 
is greater than that of the main drop. Notice that there is considerable scatter in the available data in the figure. The 
computed results from the axisymmetric LBE model are presented for wavenumbers greater than or equal to 0.47. 
Ignoring the two experimental data points of Lafrance (1975) for the satellite drop sizes that deviate considerably 
from the others, we find that the axisymmetric model is able to reproduce the drop size distribution quantitatively 
within 12%. 

The axisymmetric model has been employed to study head-on collisions of drops of radii Ri and i?2 approaching 
each other with a relative velocity U . The dynamics and outcome of colliding drops is characterized mainly by the 
Weber number. We defined by We = pi{Ri + R2)U^/ct [s^. Additional parameters that may have an influence are 
the Ohnesorge number. Oh, defined by Oh = 16/ii/V piRicr and ratios of liquid and gas densities(r) and dynamic 
viscosities (A). According to experiments ^39], it is expected that lower We collisions lead to coalescence while higher 
We to separation by reflexive action. Figures ^1 and ^] present drop configurations at We — 20 and We = 100 
respectively. Notice that at We = 20, the drops coalesce, while at We — 100, they eventually separate with the 
formation of a satellite droplet, which are consistent with experimental observations. Also notice that for the latter 
case, the temporarily coalesced drop undergoes various stages of deformation which are consistent with a recent 
theoretical analysis |40| . Additional details of these and other studies of drop collisions are given in Ref. . 
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FIG. 9: Drop sizes resulting from Rayleigh breakup of liquid cylindrical column as a function of wave number k* . Quantities 
are dimensionless. 



V. SUMMARY 



In this paper, a LB model for axisymmetric multiphase flows is developed. The axisymmetric model is developed 
by adding source terms to the standard Cartesian BGK LBE. The source terms, which arc temporally and spatially 
dependent, represent the axisymmetric contributions of the order parameter, which distinguish the different phases, 
as well as inertial, viscous and surface tension forces. Consistency of the model in achieving the desired axisymmetric 
flow multiphase behavior is established through the Chapman-Enskog multiscale analysis. The analysis shows that 
the axisymmetric macroscopic conservation equations are recovered in the continuum limit. An axisymmetric model 
with reduced compressibility effects is then developed to improve its computational stability. In this version, a 
transformation is introduced to the distribution function in the LBE such that it reduces the compressibility effects. 
Comparisons of computed axisymmetric eqiiilibrium drop formation and oscillations, Rayleigh capillary instability, 
breakup and formation of satellite drops liquid cylindrical liquid columns and the outcomes of head-on drop collisions 
with available data show satisfactory agreement. The maximum error for the frequency of drop oscillations is less 
than 5.1% and that for drop sizes as a result of Rayleigh breakup is 12%. 
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FIG. 10: Colliding drops at different times T; We — 20, Oh — 0.589, r = 4, A = 1. Time is normalized by the relative velocity 
between the drops and their diameter. Axes are in lattice units. 



[9: 

[10 

[11 

[12; 
[is; 

[14 

[is: 
[16: 
[17: 
[is: 

[19 
[20 
[21 
[22' 
[23^ 
[24' 
[25^ 
[26' 
[27' 
[28' 
[29 



Swift, W. Obsborn, and J. Yeomans, Phys. Rev. Lett. 75, 830 (1995). 
Swift, S. Orlandini, W. Obsborn, and J. Yeomans, Phys. Rev. E 54, 5041 (1996). 
He, X. Shan, and G. Doolen, Phys. Rev. E 57, R13 (1998). 
He, S. Chen, and R. Zhang, J. Comp. Phys. 152, 642 (1999). 
He and G. Doolen, J. Stat. Phys. 107, 112 (2002). 
Sussman and P. Smereka, J. Fluid. Mech. 341, 269 (1996). 
He, R. Zhang, and G. Doolen, Phys. Fluids 11, 1143 (1999). 
Inamuro, R. Tomita, and F. Ogino, Int. J. Mod. Phys. B 17, 21 (2003). 
Premnath and J. Abraham, Int. J. Mod. Phys. C. To appear (2004). 
I. Halliday, L. Hammond, C. Care, K. Good, and A. Stevens, Phys. Rev. E 64, 011208 (2001). 
B. Nadiga and S. Zaleski, Eur. J. Mech. B: Fluids 15, 885 (1996). 
Zou and X. He, Phys. Rev. E 59, 1253 (1999). 

Rowlinson and B. Widom, Molecular Theory of Capillarity (Clarendon Press, Oxford, 1982). 

Chapman and T. Cowling, Mathematical Theory of Non- Uniform Gases (Cambridge University Press, London, 1964). 
Evans, Adv. Phys. 28, 143 (1979). 

Carnahan and K. Starling, J. Chem. Phys. 51, 635 (1969). 
Qian, D. d'Humieres, and P. Lallemand, Europhys. Lett. 17, 479 (1992). 
Bhatnagar, E. Gross, and M. Krook, Phys. Rev. 94, 511 (1954). 
He and L.-S. Luo, Phys. Rev. E 55, R63333 (1997). 
He, S. Chen, and G. Doolen, J. Comp. Phys. 146, 282 (1998). 
-S. Luo, Phys. Rev. E 62, 4982 (2000). 

Guo, C. Zheng, and B. Shi, Phys. Rev. E 65, 046308 (2002). 



M. 
M. 
X. 
X. 
X. 
M. 
X. 
T. 
K. 



Q. 
J. 
S. 
R. 
N. 
Y. 
P. 
X. 
X. 
L. 
Z. 



17 



180 




180 



180 



180 



100 200 300 400 500 600 
XAxis 

{i)T = 




100 200 300 400 500 600 
XAxis 
(iii)T = 2.816 




100 200 300 400 500 600 
XAxis 
(v)T = 11.360 




100 200 300 400 500 600 
XAxis 

(vii)T= 19.904 




180 



180 



180 



S, 90 
>- 



100 200 300 400 500 600 
XAxis 

(ii)T=1.551 




100 200 300 400 500 600 
XAxis 
(iv)T = 8.879 




100 200 300 400 500 600 
XAxis 
(vi)T= 17.056 

















~^ ^ 









100 


200 300 400 
XAxis 

(viii)T = 26.549 


500 


600 



FIG. 11: Colliding drops at different times T; We = 100, Oh — 0.589, r = 4, A = 1. Time is normalized by the relative 
velocity between the drops and their diameter. Axes are in lattice units. 
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